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ABSTRACT 

Tidal streams provide a powerful tool by means of which the matter distribution of 
the dark matter halos of their host galaxies can be studied. However, the analysis is 
not straightforward because streams do not delineate orbits, and for most streams, 
especially those in external galaxies, kinematic information is absent. We present a 
method wherein streams are fit with simple corrections made to possible orbits of 
the progenitor, using a Bayesian technique known as Parallel Tempering to efficiently 
explore the parameter space. We show that it is possible to constrain the shape of the 
host halo potential or its density distribution using only the projection of tidal streams 
on the sky, if the host halo is considered to be axisymmetric. By adding kinematic 
data or the circular velocity curve of the host to the fitting data, we are able to recover 
other parameters of the matter distribution such as its mass and profile. We test our 
method on several simulated low mass stellar streams and also explore the cases for 
which additional data are required. 
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1 INTRODUCTION 

The global features of the dark matter halos of spiral galax- 
ies such as their mass, shape and extent hold key clues to 
understanding the nature of dark matter particles as well 
as galactic evolution and morphology. Simulations of cold 



dark matter tend to produce halos that are oblate (Dubin- 



ski 1994 Debattista et al. 20081 or triaxial (Frenk et al 



1988| pubinski fc Carlberg||1991| |Bett et al.||2007| > whereas 



hot dark matter models yield more spherical halos (Mayer 



et al. 2002 Bullock 2002) and baryonic dark matter such 



as cold molecular gas form disk-like halos (Pfenniger et al 



19941. It is crucial to test these theoretical predictions with 



observations and to this end, several techniques have been 
developed to probe the vertical distribution of dark matter 
around galaxies. These range from methods that are very 
local, such as measuring the density of dark matter in the 
Solar neighbourhood which is sensitive to the halo flatten- 
ing (Kuijken & Gilmorc 19891, to methods based on grav 



itational lensing (Hoekstra et al. 2004 Mandelbaum et al. 
2006). Yet, a conclusive result remains elusive. This is in 



part due to the small sample of galaxies for which these 
measurements have been made. The other cause for the dis- 
agreement in results may be the systematic differences in 
the methods themselves. For instance, Oiling & Mcrrificld 
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(2000) measured the flattening of halos based on its effect 
on the thickness of the H I gas layer in disk galaxies and find 
that their method tends to yield flatter halos as opposed to 
other techniques such as those based on warped gas disks. 
A recent study of halo shapes based on H I flaring can be 
found in |Q'Brien et al.| ( |2010a|b|c|d[ ). See |Sackett| (|1999|> for 
a review of earlier work on the subject, and |Merrifield| ( |2004[) 
for measurements of the shape of the Milky Way halo. |Alh| 
good et al. ( 2006 1 contains a more recent comparative study 



of different methods to determine halo shapes. 

Orbits in a potential are simple, yet powerful tracers 
of the matter distribution. Polar rings in early-type galax- 



ies are a good example of these (Sackett et al. 1994 Iodice 



|et al.|[2003} |Combes||2006) . In this paper, we use a simi- 
lar tracer, namely the stellar streams that are formed by 
the tidal disruption of globular clusters or dwarf spheroidal 
galaxies as they fall into the halos of larger spiral galaxies. 
The stars that make up the stream are ejected from the 
in-falling satellite due to tidal disruption and shocks dur- 
ing pericenter passages. The stream stars occupy (approxi- 
mately) two different orbits that are slightly offset from that 
of the progenitor, one that is ahead of it and one lagging be- 
hind, making up the leading and trailing arms of the stream 
respectively. The positions of the stream stars can be calcu- 
lated fairly easily from the orbit of the progenitor (see |5|). 
As the geometry of the stream depends on the orbit of the 
progenitor, which in turn depends on the potential of the 
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host galaxy, it can be used to constrain the parameters of 
the halo potential or its density distribution. Tidal streams 
in the halo have the added advantage of being far enough 
from the bright galactic components that the fitting is un- 
affected by small errors in the disk and bulge models. Tidal 
tails formed by interacting disk galaxies in major mergers 



can also be used to probe the parent halos ( Dubinski et al. 
1999| |Springel fc White||l999") , but in this study, we focus 



solely on stellar streams formed in minor mergers. 

Several stellar streams have been observed in the Milky 
Way; these include the Sagittarius (|Ibata et al.|2001b Ma- 
jewski et aL"1|2003| ), the Orphan ( |Belokurov et al.||2007| 



Grillmair|2006 1 and the Monoceros tidal streams (Newberg 
et al.|2002||Conn et al.|2005||2007[|2008[ ). Stream-like struc 



tures have also been observed around many nearby galax- 
ies, for example M31 ([Ibata et al.|[200lal [2004)), NGC 5907 



( Martmez-Delgado et al.||2008[), NGC 891 (|Mouhcine et a! 
|2010[ ) and others; refer |Martmez-Delgado et al.| ( |2010| ) for a 
recent systematic survey of streams in nearby spiral galax- 
ies. Many more remain to be uncovered by future surveys 
that are increasingly sensitive to low surface brightness ob- 
jects, enabling an extensive application of techniques based 
on them. These could provide a large sample of measure- 
ments of halo shapes and profiles, which is required if we 
are going to be able to uncover the generic properties of 
dark halos and study the possible correlations with their 
formation histories. 

There has been an ongoing effort in recent years to glean 
information from these streams about their host halo distri- 
bution. For example, the fact that the stars of the Sagittarius 
tidal stream lie on a narrow great circle on the sky (implying 
little precession) and that the stream has not been dispersed, 
initially led to the conclusion that the Milky Way has a 



spherical halo (Ibata et al. 


2001b 


Majewski et al. 20031. 


Mayer et al.| (|2002|), |Helmi 


| (|2004 


1 and IMartmez-Delgado 



et al. ( 2002 I suggested that this need not necessarily be 
the case as the debris is dynamically young and may not 
be sensitive to the shape of the halo. |Law et ak] ( |2009[ > have 
shown that the Sagittarius tidal stream is best reproduced 
in a triaxial halo and more recently, a more mildly triax- 
ial halo ( |Law fc Majewski]|2010[ ). Attempts have also been 
made to constrain the Galactic potential by fitting thin, cold 



streams found in the Milky Way halo, GDI (Grillmair & 
|Dionatos||2006|) and Palomar 5 flOdenkirchen et al.||2003l, 



and Palomar 5 (Odenkirchen et al 
with orbits integrated in a triaxial potential, but the data 
were found to be insufficient to discriminate between pos- 
sible solutions, although triaxial models seem to be favored 
over spherical ones (Lux 2010 1. However, by fitting the GD-1 



stream with orbits in a logarithmic, axisymmetric potential, 



Koposov et al. (20101 have estimated the flattening of the 
Milky Way potential at galactocentric radii near R ~ 15kpc 
at 0.9 and a lower limit for the flattening of its halo potential 
at 0.89. 

In general, there are two main challenges that have so 
far hindered the successful application of tidal streams in 
constraining halo shapes. One is that tidal streams do not 
exactly delineate individual orbits in the galactic potential 
| |Eyre fc Bi rmey 2009a b| and treating them as such may 
result in incorrect estimates. We overcome this hurdle by 
fitting a given stream with corrected sets of points com- 
puted from the progenitor orbit. The other limitation is that 
in most systems, only projected coordinates of the stream 



on the sky are available. In closer systems such as the An- 
dromeda galaxy (M31), line of sight velocities are also avail- 
able and distances to globular clusters or dwarf galaxies that 
may be the stream progenitors are measurable using the tip 



of the Red Giant Branch (TRGB method - see McConnachie 
|et~aTl2004l In previous studies using tidal streams, the ap- 
proach has been to try to reproduce the observed streams us- 
ing a few N-body simulations (typically less than 10). It has 
been found that the phase space information is insufficient 
in reasonably constraining the halo parameters using this 
approach. The N-body technique is severely handicapped 
by the fact that it is only possible to explore a tiny frac- 
tion of the full parameter space with such simulations. We 
redress this shortcoming by adopting a statistical approach 
sampling the parameter space of possible progenitor orbits 
and halo parameters using a Markov chain Monte Carlo rou- 
tine. The sampling yields distributions of the halo structural 
parameters that peak at their most likely values. We find in 
our tests on simulated streams that it is possible to estimate 
halo shapes using this method even with limited phase space 
information. We also find that it is possible to recover line 
of sight distances along streams in cases where they are not 
available. In this paper, we explore how much and what in- 
formation is required to uniquely estimate the orbital and 
potential parameters using this statistical approach and how 
the estimates improve with additional data. Later contribu- 
tions in this series will use the technique developed here to 
actual observed systems. 



2 METHODOLOGY 

Consider a galaxy whose disk lies in the XY plane. We as- 
sume that the mid-plane of its bulge and halo coincide with 
this plane. We also assume that the halos we study are static 
and have an axisymmetric density profile, even though the 
density profile of the dark halo is generally considered to be 



a triaxial ellipsoid ( Dubinski fc Carlberg||1991 Jing & Suto 
|2002||Bailin fc Steinmetz|2005[ ). The flatness of the distribu- 
tion is given by the ratio of the polar to the equatorial axis 
(c/a) and the ovalness by the equatorial axis ratio (b/a). In 
the most general models, which involve a superposition of 
ellipsoids, c/a and b/a vary with radius. 

However, for the present contribution, we make the 
simplifying assumption that the equatorial axes are equal 
(b/a — 1), and that the halo has a flattening of q — c/a (we 
discuss both flattening in the potential and density distribu- 
tion). Studies have shown that this is a reasonable assump- 
tion for the halos of spiral galaxies, especially far from the 



galactic disk where it becomes spheroidal (Dcbattista et al. 
2008} |Kazantzidis et al.||2010[ ) 



Throughout the paper, we analyse the specific case of 
a host galaxy that is viewed edge-on, although the method 
can be very simply adapted to systems that do not have that 
geometry. What we have access to are the projected coordi- 
nates (x,z) of a tidal stream in its potential, and the width 
of the stream providing an uncertainty of a on these coor- 
dinate values. The tidal stream has two tails: the leading 
and the trailing tails. In reality, each tail consists of several 
different yet very similar orbits, the stars along these having 
similar values of energy and angular momentum. The slight 
difference in these orbits show up in places in the form of 
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bifurcations and small protrusions from the main stellar fea- 
ture. In this paper, we neglect these tiny features and only 
consider the longest, contiguous structures that we observe. 
The method is only applicable in its current simple form 
to streams of low mass satellites ( 10 8 Mg), for which the 
self-gravity of the stream is negligible. This allows us to ig- 
nore dynamical friction which would have to be taken into 
consideration for more massive and heavily disrupted sys- 
tems. 

The method for fitting the streams is briefly as follows. 
Different values of the halo potential parameters and the or- 
bital parameters of the progenitor are tested, each set of val- 
ues corresponding to a point in the multidimensional param- 
eter space. For each set of parameters considered, an orbit 
is integrated and a set of points is calculated that represents 
the stream that would be formed by a satellite on the orbit 
(see [JHJ; we refer to this as the trial stream corresponding to 
the point in parameter space under consideration. Compar- 
ing the trial stream to the observed stream (or the N-body 
generated test stream in the present work), we measure how 
likely the parameters are to be the ones that generated the 
observed stream. We sample different values of the parame- 
ters using a Markov chain Monte Carlo (MCMC) algorithm, 
in order to deduce the distribution of parameter values that 
are consistent with the data. For the galactic halo potential, 
we use two models: a purely logarithmic halo for its simplic- 
ity and a more realistic multiple-slope power-law model. 

In Q we explore how the quality of the estimate de- 
pends on the quantity and type of information available. 
We discuss cases for which the streams are shorter or have 
fewer turning points. We also consider the different kinds of 
data that may be available and their effect on the accuracy 
of the estimate. For instance, line of sight velocities and/or 
distances at some points along the stream or the rotational 
velocity curve are available for nearby systems. The method 
is tested on these various cases with pseudo-data generated 
by N-body simulations. 



3 THE FITTING ROUTINE 

In this section, we first briefly discuss the general principles 
of maximum likelihood estimation and parallel tempering, 
and then describe how we apply it to fitting tidal streams. 

Maximum likelihood parameter estimation is a robust 
method of determining the parameters of a model that max- 
imize the likelihood of a data set. The likelihood of a set of 
parameters is the probability of obtaining the given data for 
those parameter values. Suppose the data set consists of N 
independent observations of a random variable x which has 
a probability distribution function f(x;<d) that depends on 
a set of k parameters 6 which are to be estimated, then the 
likelihood is calculated by the product of the values of the 
distribution function for each data point ( Gregory |2005 1: 



C(xi,x 2 ,...,x N \@i,e 2 ,...,® k ) =JJ/(a:<; 61,62, ...,0fc).(l) 

i=l 

The estimated values of the parameters 6 are the ones for 
which this likelihood is maximum. A basic and simple algo- 
rithm to find this set of parameters is known as the Metropo- 
lis Algorithm, which probes the parameter space in an effi- 



cient way to locate the parameters for which the likelihoods 
are high. The sampling is done by a chain that walks through 
the parameter space. 

For a single chain, the walk through the parameter 
space basically involves choosing a random point that we 
consider to move to, and then moving to it or staying put 
based on a likelihood criterion. Such a chain is known as a 
Markov chain Monte Carlo (MCMC). A chain is set off at an 
initial random point Qt=o in the parameter space (t being 
the time step of the MCMC), and we calculate the likeli- 
hood at this point, say £(6t). Then we consider another 
point 6' (the trial point) in the parameter space, from a 
proposal distribution and calculate its likelihood £■(&'). In 
our algorithm, the proposal distribution is a normal distri- 
bution centered on the current point. If the trial point falls 
beyond the plausible range of the parameters, then we set its 
likelihood to be arbitrarily low. We calculate the Metropolis 
ratio, r = £(6t)/£(6'). If this is greater than 1, then we 
take the trial point as the next point in the chain. If r < 1, 
then we pick a uniform random number U between and 
1. If U ^ r, then we take the trial point as the next point 
(6t+i = 6'), else we stay at the same point (6t+i = 0t). 
Proceeding in a similar fashion, we explore the full parame- 
ter space. 

3.0.1 Parallel Tempering 

In many cases, especially in multidimensional parameter 
space, there are usually many local likelihood maxima, 
whereas we seek the global maximum. Using only a single 
chain, the algorithm often becomes stuck in a local maxi- 
mum. To avoid this, multiple chains of different "temper- 
atures" are used to step through the parameter space. Ef- 
fectively, this means that for a chain of temperature T, its 
likelihood is raised to the power of the inverse of its temper- 
ature, i.e.: 



(2) 



where C\ corresponds to the likelihood at a temperature of 1 
(given by equation [T]in general, which reduces to equation|4] 
for the present study) and /3; = 1/Tj. T, is the temperature 
of the i th chain. As T ranges from 1 to infinity, /3 takes values 
from 1 to 0. The higher temperature chains are less sensitive 
to the likelihood, and take larger strides through the param- 
eter space. Suppose we use n parallel chains and swap their 
states after every n a steps on average. The higher temper- 
ature chains pull out the colder chains from local maxima 
they may be stuck in. The algorithm for the swapping of 
chains is as follows: 

(i) At each step, choose a uniform random number Ui 
between and 1. A swap is proposed if U\ ^ l/n s . 

(ii) If a swap is proposed, we consider swapping the states 
of chain i and i + 1, where i is a random integer between 1 
and n — 1. 

(iii) Calculate the swapping probability, 

r ft (e t ,, +1 )£ ft+1 (e M ) 



£ft(e M )£/3 1+1 (e t , l+ i) 



(3) 



where Q t ,i and 6t,i+i are the current points on the i* 



and (i - 



chains respectively, and £^(6^4+1) is the like- 



lihood of 6t,j+i if it were on the i chain. 
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(iv) Choose another uniform random number U2 between 
and 1. Accept the swap if U2 < r. 

(v) The maximum corresponding to the coldest chain T = 
P = 1 is the global maximum. 

The efficiency of the algorithm is sensitive to the size of the 
proposal distribution. If it is too small, then most of the 
trial points are accepted and the MCMC is slow to sample 
the full parameter space. If it is too large, most of the trial 
points are rejected and though the MCMC may make large 
jumps in parameter space, it could become stuck at a cer- 
tain point for long. For the present work, optimal sizes of the 
proposal distributions for the parameters were found exper- 
imentally, by requiring that the acceptance rate lie between 
15% and 40%. Note that it is possible, however, to automate 
the process such that the routine finds an optimal size for 
the proposal distribution (Gregory 2005). It is important to 
check that the MCMC converges on a solution. One way of 
checking for convergence is to sample the parameter space 
with chains that start at different initial points and to check 
if they yield similar solutions. Once convergence has been 
achieved (i.e. a chain is well-mixed), the distribution is in- 
dependent of the number of MCMC steps and is said to be 
stationary. The initial steps of the chain are discarded to for- 
get the starting point. This is known as burn-in and is typ- 
ically around 10, 000 steps in our test cases. The remaining 
steps of the chain form a sample distribution of the parame- 
ters. The estimated value of a parameter corresponds to the 
peak of its marginalized distribution. The accuracy of the 
estimate is equal to the standard deviation of the distribu- 
tion, if it turns out to approximate a Gaussian distribution. 
However, the parameter distributions may turn out to be 
multimodal or otherwise complex. As the specific details of 
the parameter distributions vary from stream to stream, we 
defer the discussion of accuracies and convergence tests to 
when we apply the method to real, observed streams. For 
a detailed discussion of Parallel Tempering and maximum 
likelihood estimation, see |Gregory|prjQ5| . 



3.1 Estimation of Halo Parameters 

By defining the parameter space and the likelihood calcula- 
tion for the parameters, we can use the above method to find 
an orbit that best reproduces a given stream (after the re- 
quired corrections explained in ^5|, thereby also estimating 
the model parameters that we wish to determine. A point 
in the parameter space corresponds to an orbit, which is 
parametrized by its initial position and velocity components 
(xo, yo, 20, v x<3 , v yo , v ZQ at t = 0), in addition to the poten- 
tial or density distribution parameters. Depending on the 
information available, few of these initial phase space co- 
ordinates are known. The remaining orbital parameters are 
set to be free parameters with physically reasonable ranges. 
The best fit orbit corresponds to the one with the highest 
likelihood. The main processes involved in the algorithm are: 
(1) Stepping through the parameter space, (2) Integrating 
orbits corresponding to the points in the parameter space, 
(3) Generating a corrected set of points (trial stream) cor- 
responding to an orbit and (4) Calculating the likelihood of 
these trial streams. 

As the several MCMC chains step through the param- 
eter space, we integrate orbits at each step with a fourth- 



order Runge-Kutta integrator starting from an initial point 
defined by the parameters. For fitting the streams, we cal- 
culate a set of n points {S(n)} via a correction mechanism 
(detailed in i|5| from the integrated orbit. Assuming a Gaus- 
sian distribution for the parameters, the likelihood of a trial 
orbit is calculated as: 



2tt 



(4) 



where TV is the number of data points and cr, is the uncer- 
tainty in each data point, which includes both the width 
of the stream as well as measurement errors. Di is the dis- 
tance of each data point on the observed stream from the 
corrected set of points, i.e. the distance of a data point to 
its closest point on the trial stream. As not all of the phase 
space information is available, the distance Di only includes 
the distances along the coordinates for which data are avail- 
able. For instance, if the projected positions Xi,Zi and line 
of sight (l.o. s) velocities v Vi are observed, then the exponent 
in equation [4] reduces to: 



2a? 



E 



(Xi - x c ) 2 + (Zi — z c ) 2 (v 



lol 



+ 



2ai 



(5) 



where x c , z c , Vy c are the closest points on the trial stream to 
the data point Xi,Zi,v yi and a XZi ,a Vi are the uncertainties 
in position and velocity. The <Tj in the normalizing coefficient 
of the exponent in equation [4] is the product of the uncer- 
tainties in the different phase space coordinates, which in 
this example reduces to a XZi a Vi . 

In our tests, we used four parallel chains to search the 
parameter space, the inverse of the temperature of the i th 
chain being /3, = 1/i 2 . The states of the chains are con- 
sidered for swapping after every thirty steps on average. 
At the beginning of the parallel tempering, all the chains 
are initialized by setting the free parameters to physically 
plausible random values. We checked that the results were 
insensitive to these starting values. Given the stream data, 
the likelihood of the orbit corresponding to the initial point 
is calculated with equation [4] This is the likelihood of the 
starting point for the coldest chain, i = 1. The likelihoods 
for the same point on the hotter chains are calculated using 
equation [2] 

Having calculated the likelihoods for the initial point, 
the next point on each MCMC has to be chosen. For this, 
a trial point is chosen for each chain. This can be done in 
many ways. In our case, if Ot is the current value of the 
parameter O, then the next point for consideration 0' is 
drawn randomly from a normal distribution centered on t . 
The width of the proposal distribution, i.e. the size of the 
step, depends on the expected range of the parameters. It 
also helps to use different step sizes for the different chains 
to enable sampling the parameter space at various scales. 
The algorithm is run until the chains are well mixed, which 
could take from approximately one hundred thousand to a 
million steps (depending on the stream and halo model) , and 
the marginalized distributions of the parameters are drawn 
from the coldest chain. 
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4 TESTING WITH ORBITS (NOT STREAMS) 

To get a heuristic idea of the effectiveness of the method, 
we first investigated how well it constrains potential param- 
eters using ideal orbits in both logarithmic halos and double 
power law density halos. The primary aim was to explore if 
and when the projected positions of an orbit are sufficient 
to estimate the density or potential flattening and how ad- 
ditional information would affect the accuracy of the esti- 
mation. 



4.1 Orbits in a Logarithmic Potential 

Consider the set of orbits shown in Fig ure [Tj These were in- 
tegrated in a logarithmic halo given by ( Bir mey fc Tremaine| 
20(1]): 



1 , 

$haio = -V ln 



-R? + R 2 H — .. 



(6) 



in cylindrical coordinates, where q$ is the potential flatten- 
ing, Vb is the circular velocity, and R c is the core radius. For 
this initial test, we neglect the contribution of the bulge and 
the disk to the potential. The orbital and potential parame- 
ters of each of these are listed in Table[l] A logarithmic halo 
has the advantage of having only three parameters, thereby 
minimizing possible degeneracies between the different mod- 
els. Using only the projected positions of these orbits, we find 
that the flattening q,p and initial line of sight distance of the 
orbit from the center of the galaxy yo are easily constrained 
as shown in Figures [2] and [3] whereas the circular velocity is 
degenerate and cannot be constrained with only positional 
information (Figure [4J. It is to be noted that there may be 
a sign discrepancy in the distance estimate, as orbits with 
either value of yo will be identical in projection. Therefore, 
we restrict the fitting routine to positive yo space and only 
the magnitude of the progenitor distance is recovered. The 
orbits in this set are fairly long and have more than two 
turning points. However, several observed streams are much 
shorter. In order to analyze the effect of the length of an 
orbit and more importantly, the number of turning points 
of the orbit, on the estimate of the flattening, we fit shorter 
versions of the orbit AS1 (Figure[5|. Not surprisingly, having 
fewer turning points on an orbit causes a spread in the es- 
timated value of q$. Adding more information to the fitting 
(line of sight velocities for BS1 and line of sight velocities 
and distances y for CS1) results in much more accurate es- 
timations of q^ as shown in the bottom panels of Figure [5] 
As for the distance to the progenitor, yo, it is strongly con- 
strained (with a sign discrepancy) for BS1 even only with 
the projected positions and for CS1, if kinematic informa- 
tion is provided (but not with only the projection), making 
it the most easily constrainable parameter. 



4.2 Orbits in a Spheroidal Halo 

A more general halo is described by a double power law 
density distribution ( Dehnen fc Binney]|1998 \ 

where s is the ellipsoidal coordinate 
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Figure 1. Projection in the XZ plane of orbits integrated in a 
logarithmic potential using a Runge-Kutta scheme. The dots indi- 
cate the positions on the orbits that are used as data points in the 
fitting. The top, middle and bottom panels show orbits in spheri- 
cal, oblate and prolate potentials respectively. The parameters of 
each orbit are listed in Tabled 
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Figure 2. Estimation of <j0 for the projected orbits shown in 
Figurejl] The input value of in each case is shown. This distri- 
bution is drawn from 100,000 steps of the coldest MCMC chain. 
The excellent correspondence with the input values shows that 
the shape of a logarithmic potential can be accurately recovered 
from the spatial projection of some orbits within it. 
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Table 1. Series of Logarithmic Orbits: The following are the potential and orbital parameters of the orbits shown in Figure [T] q^^Rc 
and Vq are the potential flattening, core radius and circular velocity respectively. xo,yo,zo a- re the initial positions and v xo ,v yo ,v zo the 
initial velocities used. 
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Figure 3. Estimation of yo f° r the projected orbits shown in 
Figure^ The input value of yo in each case is shown. This distri- 
bution is drawn from 100,000 steps of the coldest MCMC chain, 
and clearly peaks near the input value. The fitting routine only 
considers positive values for yo as orbits that are symmetrical 
about the XZ plane will have identical positions (so there may 
be a sign discrepancy in the estimated value of yo). 



Figure 4. Estimation of Vo for the orbits shown in Figure [l] The 
input value of Vo in each case is shown. This distribution is drawn 
from 100,000 steps of the coldest MCMC chain. We see that the 
peaks of the distribution vary greatly from the input value, which 
implies that the orbits are not uniquely dependent on the circular 
velocity. 



s = [R + z /q p 



2 , 2\l/2 



(8) 



The model parameters that are unknown are the central 
density po, the inner and outer slopes 7, /3, the scale radius 
ro and the flattening in density q p . We calculate the poten- 
tial due to this distribution by multipole expansion using 
an algorithm similar to 'GalPot' (De hnen fc Binney||1998| 
Binney fc Tremaine|2008 |. The truncation radius rt is fixed 
at 1000 kpc, but its precise value does not affect orbits at 
the radial distances under consideration. We find that the 
inner slope 7 cannot be constrained by orbits that are far 
away from the center of the galaxy and that small variations 
in 7 do not affect the fits. Hence, we adopt a fixed value of 
7=1, equivalent to the central power-law slope of a NFW 



profile ( |Navarro et al.||l996 l, for the orbital fitting and set 
the remaining parameters to be free. Figure [6] shows the dis- 
tributions in the flattening q p obtained by fitting the projec- 
tions of three orbits in this power-law spheroidal halo. Even 
though the distributions peak at the right values of q p , they 
are much more spread out than the previous distributions in 
q$. This is to be expected, as orbits respond directly to the 
potential, and these are rounder than their corresponding 
mass distributions. For a logarithmic potential, at distances 
much larger than the core radius, the relation between q^ 
and q p can be approximated as (Binney & Tremainc 2008): 



1 - Q4> w 3 (1 - Qp) 



(9) 
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Figure 5. Fitting shorter versions of the orbit AS1 (q^, = 1). 
BS1 has two turning points and CS1 only one. The middle panels 
show the <j0 distributions obtained by fitting only the projected 
positions of these orbits. The spread in the distribution increases 
with decreasing turning points. The bottom panels show the q^ 
distribution obtained by adding more information to the fitting: 
the line of sight velocities v y for BS1, the distances y and line of 
sight velocities v y for CS1. 



The fitting mechanism also turns out to be extremely useful 
in approximating the line of sight distances along the orbit 
as revealed by a grayscale plot of the distances along the 
various trial orbits on the coldest MCMC chain (Figure f7|. 



5 TESTING WITH STREAMS IN A 
SPHEROIDAL POTENTIAL 

Having demonstrated the power of the technique in con- 
straining the parameters of a density distribution by fitting 
only the projected positions of orbits, we test its ability to 
robustly estimate the same by fitting streams which, as men- 
tioned in [jl] do not delineate the orbit of the progenitor or 
any other exact orbit in the potential. A consistent and fast 
mechanism is required to derive the positions of stream stars 
from the progenitor's orbit, without using N-body integra- 
tion. The stars which make up the stream are those which 
were tidally ripped from the satellite during pericenter pas- 
sages. Based on this, it is possible to formulate a simple 
correction mechanism that maps a given progenitor orbit to 
the coordinates of its tidal tails. 

Consider the stream shown in Figure |SJ whose projec- 
tion in the V X V Z plane is shown in the right panel of Figure 
[9] (it is the same as stream B in Figure |10| which was in- 
tegrated in an oblate halo). The grey dotted curve is the 
orbit of the progenitor, the remnant of which is the concen- 
trated sphere. The progenitor's orbit corresponds to a set 
of orbital parameters (the progenitor's current position and 
velocity components) and potential parameters. The stars 
that lie on the tidal tails escaped the progenitor at earlier 



times, from two regions of the satellite: one that is closest to 
the center of the host galaxy and one that is farthest from 
it. We refer to these points as the inner and outer escape 
points respectively, which approximate the inner and outer 
Lagrange points in a restricted three-body problem. There- 
fore, the trajectory of a stream star can be integrated with 
its initial position offset from a certain point in the progeni- 
tor's past orbit by a certain distance r cuto ff, offset outwards 
for trailing tail stars and inwards for leading tail stars. This 
is shown diagrammatically in Figure [8] Starting at the cur- 
rent position and velocity of the progenitor P, its orbit is 
integrated backwards in time, marked by the grey dots. For 
any point on its backward orbit, Q, at time tQ, the inner 
and outer escape points are approximated as tq — r cu toff 
and tq +r cu toff, indicated respectively by the blue and red 
dots at Q. These provide the initial positions for the orbits 
of stars that escape the satellite at Q. For their initial veloc- 
ity components, we use the velocity of the progenitor orbit 
at Q. With these initial phase space coordinates, we inte- 
grate forward for the same amount of time £q . These orbits 
are indicated by the blue and red dotted curves, with their 
final points on the leading and trailing tails at A and B re- 
spectively. Repeating this process for several points on the 
backward integrated orbit (say, for every 50 Myr) yields a 
set of points which lie closely on the tidal tails of the stream 
(red and blue dots in Figure |9| , which we shall refer to as 
the corrected points for a given set of parameters. 

We find empirically that the offset radius for a point 
Q at a radial distance r can be calculated approximately as 
2.88 times the theoretical Jacobi radius: 



1/3 



rcutoff = 2.88 x — | 

\3M(r), 



(10) 



where m sat is the mass of the satellite and M(r) is calcu- 
lated from the circular velocity V c of the host galaxy at r as 
M(r) = rVc(r)/G, where G is the gravitational constant. 
For a perfectly spherical potential, M(r) is equivalent to 
the mass contained within radius r. For non-spherical po- 
tentials, it is only a crude approximation to the total mass 
inside r, but as we show below, the correction we obtain us- 
ing Tcutof f calculated in this manner is sufficiently accurate 
for our purposes. 

We find that by fitting the stream data with these cor- 
rected points it is indeed possible to recover the parameters 
of the potential as well as the orbit. It is interesting to note 
that through this correction mechanism, one has more infor- 
mation on the progenitor's past orbit than one would have 
with only the local orbit of the progenitor (i.e. the red and 
blue dotted curves in Figure [9] contain more information 
than the grey curve does). In this sense, the tidal stream 
retains some memory of its infall. It is important to note 
that although we do not correct for velocities at the escape 
points, the velocities of the corrected points also lie along 
the velocities of the tidal tails as seen in Figure [9] i.e. the 
correction is valid in all of phase space. This is true of all the 
streams we tested, possibly because we use low mass satel- 
lites. If a similar correction mechanism is used for extending 
the present technique to more massive streams, we expect 
that a corresponding correction in velocity space would be 
required. 

The backward integration time depends on the age of 
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Figure 6. Estimation of q p for the orbits shown in the top panel. The host galaxy models in these examples are a one-component 
ellipsoidal halo (Equation[7jl. The input value of q p in each case is marked in red. These distributions are drawn from 2,000,000 steps of 
the coldest MCMC chain. 
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Figure 7. Recovery of the three dimensional structure of orbits with only their projections on the XZ plane. The grayscale image shows 
the density of the line of sight distances along the trial orbits considered in the coldest MCMC chain used in the fitting of the orbits 
shown in Figure [6] The red line marks the actual distance along the input orbits. This figure shows that coarse estimates of the distances 
along the orbit can be obtained. 



the stream and is usually unknown. It is also assumed in 
the present discussion that the progenitor is identifiable and 
its position on the sky is known, but more often than not, 
the satellite is completely disrupted or obscured and it is 
impossible to identify the progenitor (for example, the NGC 
5907 and the Monoceros streams). The time for backward 
integration can be found empirically, running different fit- 
ting routines with different backward integration times for 
the correction (it only need be accurate to around 1 Gyr). 
Similarly, several fitting routines can be tried assuming the 
unidentified progenitor to be at different positions along the 



stream. However, in all the tests described in this paper, the 
progenitor, its initial mass and backward integration time, 
were taken to be known quantities. 

Figure [10] shows nine streams on which this stream- 
fitting method was tested (note that we are now fitting 
streams, not orbits as in @. These were generated with 
gyrfalc( ON), an N-body integrator based on a fast and mo- 
mentum conserving tree-code (Dehnen 2000 20021, as im- 



plemented in the stellar dynamics toolbox NEMO (Teuben 



1995h. For the host galaxy, realistic models were used that 
comprised a spheroidal halo, a spheroidal bulge, and expo- 
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Figure 8. Correction Mechanism: Calculating a pair of corrected 
points on the leading and trailing tails corresponding to a point on 
the past orbit of the progenitor. The black dots are particles from 
an N-body simulation of the stream. The grey dotted curve shows 
the backward integrated orbit of the progenitor, starting from its 
current position at P. Q is an arbitrary point on this orbit. The 
blue and red dotted curves are the orbits of stars which escape the 
progenitor system at Q from the inner and outer escape points 
respectively. The end points of these orbits lie on the leading and 
trailing arms of the stream, marked as A and B. The velocity 
components at Q of the blue and red orbits are the same as that 
of the grey orbit (with the opposite sign). This process is repeated 
for several points along the grey curve to obtain a set of corrected 
points like A and B that demarcate the stream path (see Figure 



nential disks (thick disk, thin disk and ISM) similar to the 
ones described in Dehnen & Binney 1 1998[ ) for the Milky 
Way. Except for the flattening of the halo, the model of the 
host galaxy for each case is one of the two models described 
in Table [2] The halo flattening and the initial conditions of 
the orbits of the satellites for each stream along with its in- 
tegration time are listed in Table [3] Spherical King models 
of 100,000 particles were used for the N-body satellite sim- 
ulations. See Table [3] for the masses, tidal radii and central 
potentials of these. We have limited the present study to 
satellite masses in the range of 5.0 x 10 6 to 1.0 x 10 s M , 
as they seemed to be ideal candidates for this technique 
rather than more massive and heavily disrupted satellites, 
for which dynamical friction plays a significant role. For the 
integration with gyrfalc(ON), a Plummer softening kernel, 
along with a softening length of 0.03 kpc, and an opening 
angle tolerance parameter of 0.6 (the default) were used. 
We make extensive use of the work of Dehnen & Bin- 



ney ( 1998 ) for the calculation of the potential, especially 



for the bulge and disk components. We assume that these 
baryonic components of the system are well understood and 
their contribution to the potential is kept fixed, so that the 
only parameters to be estimated are those of the dark halo. 
As the streams are distant, the precise details of the inner 



mass distribution are not critical and the estimates are not 
affected by small variations in the baryonic component. The 
potential due to the trial density distribution of the halo is 
calculated using multipole expansion. The calculated halo 
potential and its first derivatives (w.r.t lnr and |z/r|) are 
stored on a grid in lnr and \z/r\. Adding to these, the po- 
tential and accelerations due to the bulge and the disk using 
GalPot, we have the total potential and acceleration at each 
of the grid points. The acceleration at any point is obtained 
by interpolating with a 2-dimensional, fifth-order spline. 

In the following sections, we describe the results of fit- 
ting the streams with different kinds of information. 

5.1 With projected positions alone 

In this section, we discuss how well the various parameters 
of the halo are constrained by fitting only the projection of 
streams on the sky, as this is the only information that is 
usually available for very distant systems. Figure \TT\ shows 
the distributions of the density flattening q p obtained by 
fitting the pseudo-streams with corrections to trial orbits 
described in the previous section. The number of attempted 
orbits in the fitting routines are in the range of 500, 000 
to 1,000,000, running the fitting algorithm until the chains 
are well mixed and reasonable fits of the test streams are 
found. As can be seen, the estimates fall within reasonable 
ranges of the true values of q p (indicated by the red lines), 
but the accuracy of the estimates vary from case to case. 
The differences in the accuracies arise from the differences 
in the streams themselves, as some configurations are more 
degenerate, in that many parameter choices can reproduce 
the same observations. However, even in the poor cases, the 
estimated q p still provides a useful indicator of the shape 
of the halo. For example, for stream I, made in a prolate 
halo of q p = 2.0, there are two peaks in the distribution, the 
true value of q p being at the smaller peak. The trial streams 
that are generated at both the peak values of q p resemble 
the given stream, but streams generated at the smaller peak 
(with q p = 2.0) show extra features such as small bifurca- 
tions that are similar to the ones seen in stream I, whereas 
streams generated at the larger peak (with q p ~ 1.6) do not 
show these. Thus, the streams at q p = 2.0 resemble the test 
stream in its finer details. These smaller features were not in- 
cluded in the fitting as we deemed them difficult to observe, 
especially in distant systems, but if observed, they could be 
used to obtain better estimates or distinguish between mul- 
tiple solutions. Nonetheless, the distribution in q p for stream 
I is limited to the prolate region and the corresponding dis- 
tribution in the potential flattening q^ is narrower (approxi- 
mately, q p = 1.6 = q$ = 1.2 and q p = 2.0 = q^, — 1.33, using 
equation[9|. It should be noted that for streams B and E, the 
true value of the inner slope 7 of the halo is approximately 
1.28, but they were fitted in a model with a fixed 7 = 1.0. 
This was primarily undertaken to check if the difference in 
the inner slope affects the estimation of q p , and it does not, 
confirming our assumption that these streams are insensi- 
tive to the finer details of the inner matter distribution. On 
the other hand, for fitting stream H, we set the inner power 
slope 7 to be a free parameter, with lower and upper bounds 
at -1 and 2.0. This does not seem to deteriorate the estimate 
of q p as was initially expected, so 7 can be set to be a free 
parameter in the application of the technique to real stellar 
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Figure 9. Correcting for tidal tails: The left and right panels show stream B in the XZ and V X V Z planes, respectively. The black dots 
are particles from an N-body simulation of the stream. The grey dotted curve is the orbit of the progenitor, the remnant of which is 
the concentrated sphere. The blue and red dots are the corrected points for the leading and trailing arms, respectively. The very close 
correspondence between the computationally-expensive N-body simulation and the locus of the "corrected stream" points is evident. 



streams. However, 7 itself is not found to be constrained by 
fitting the projections of these streams. 

Another significant parameter of the halo is its mass, 
which determines the circular velocities at large radii. In or- 
der to see how well our technique estimates this quantity, 
we plot the distribution of the circular velocity at 50 kpc, 
V50, shown in Figure [T2| calculated from the potential pa- 
rameters at each step of the coldest MCMC chain. As seen 
in the case for logarithmic orbits, these distributions are de- 
generate and it is not possible to constrain the halo mass 
with only positional information of a stream. We find that 
the remaining parameters of the halo, i.e. the outer slope 
P, the central density po and the scale radius ro can have 
different values for the same value of the flattening q p to 
produce streams that look identical in projection. The dis- 
tributions of these parameters are degenerate (as shown in 



measure with current instrumentation, which provides in- 
formation on the inner mass profile of the galaxy. Figure 



Figure 13 for (3) and the projection of a stream alone cannot 



constrain them. The combination of these parameters result 
in different circular velocities, making the distribution of V50 
degenerate. 

As for the orbital parameters, the line of sight distance 
yo of the progenitor from the center of the galaxy is well- 
constrained (although there may be a sign discrepancy, as 
was seen in the case of pure orbits), whereas the line of 
sight velocity or the tangential velocity components of the 
progenitor cannot be recovered. 



5.2 Adding circular velocities 

For many spiral galaxies, inner rotational velocities are avail- 
able (for instance from H I kinematics) or are feasible to 



14 



shows the effects of adding the rotational velocities (here 
assumed to extend up to 30 kpc) to the projected positions 
for stream B (to take a particular example). Comparing the 
corresponding distributions for the stream in Figures [T2| and 
T3j we see that there is a marked improvement in the esti- 
mate of the outer power slope (3. This is, as explained in 
the previous section, due to the correlation of f3 and the 
circular velocity parameter V50, the latter being estimated 
accurately with the provision of the inner circular velocities. 
The peaks of the distributions in q p in Figure [14] and Figure 
TTJare at similar positions, but the former is less smooth as 
a different proposal step size was used to fit the stream in 
this case. Though not shown here, the degeneracies in the 
scale radius ro and the central density po are also removed 
and we obtain accurate estimates of these parameters. More- 
over, distance to the progenitor and its velocities can also be 
constrained when the rotational curve is added to the fitting 
data. As seen in the distance parameter, there may be a sign 
discrepancy in the estimate of the line of sight velocity v vo 
of the progenitor, but not in the estimates of its tangential 
velocity components. 



5.3 Adding line of sight velocities 

Using stream G as an example, we illustrate the effects of 
adding kinematic data to the projected geometry of the 
stream (but without circular velocities) on the estimates of 
q p and Vso- Figure [15] shows the q p distributions in each 
case in the top panels and the V50 distributions in the bot- 
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Figure 10. The projection of test streams generated by N-body integrations using gyrfalc(ON) in a realistic multi-component galaxy 
model with bulge, disk, thick disk, ISM and halo components. The q p value indicated on each panel is the flattening of the halo used in 
each simulation. The streams in the top, middle and bottom panels were generated from satellites with initial masses of 5.0 X 10 6 M Q , 
5.0 X 10 7 Mq, and 1.0 X 10 s Mq respectively. The green dots are the points taken along the streams to fit trial streams to. 



torn panels. In (a), only the projected positions of stream 
G on the XZ plane are used in the fitting routine. In ad- 
dition to these, the line of sight velocity of the progenitor 
is provided to obtain the distributions for (b). To make the 
distributions in (c), we added more line of sight velocities 
(at five points along each tail, within 20 kpc from the pro- 
genitor). This shows that by providing velocity information 
it is possible to estimate V50 or the mass of the halo. The 
added information also improves the estimate on q p . As in 
the case where we provided the rotation curve, the progen- 
itor's velocity components and its line of sight distance are 
well constrained. 



5.4 Fitting in a logarithmic potential 

We also investigated how well the test streams in Figure [To"| 
(simulated within a full multi-component galaxy model) can 
be fit with a simple axisymmetric logarithmic halo with no 
stellar components. Figure [l6| shows the distribution in q^, 
thus obtained. The green lines indicate the potential flat- 
tening q<f, that correspond to the input density flattening q p 
values. We find that the estimates are remarkably accurate. 
This exercise suggests that a logarithmic halo can be used 
as a good approximation to the more complex input axisym- 



metric model, at least in order to constrain the flattening of 
the halo. 



5.5 Shorter streams 



We have seen in §4.1| that the estimation of parameters is 
sensitive to the number of turning points of the orbit. To 
explore this aspect of the problem for streams, we consider 
stream B at much earlier stages of its infall, where its tidal 
tails are substantially shorter. In Figure |17[ it is at 5 Gyr 
into its infall and the tidal arms form a parabolic curve, 
much like the structure observed to the north west of M31 

it is at 2 Gyr in its 



18 



( Carlberg et al.||2011[ ). In Figure 

evolution and resembles the Palomar 5 stream ([Odcnkirchcn 
et al.||200T |. For both these cases, we consider five possible 
scenarios of available information. Case A is where we only 
have the projected positions of the streams on the sky. We 
see that the distributions of q p are much more spread out 
than for stream B, the degeneracy being much more for the 
Palomar 5-like stream. For Case B, the line of sight velocities 
at certain points (red squares) are also provided for the fit. 
In Case C, we assume that no kinematic data is available, 
but the distance to the progenitor is known. It is possible to 
measure this quantity for many nearby systems as well as 
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Figure 11. Estimation of the density flattening q p with only projections of streams. Shown here are the distributions of q p obtained by 
fitting the projection of streams shown in Figure [To] The number of progenitor orbits attempted vary from 500,000 to 1,000,000 from 
stream to stream. The red lines indicate the input q p for each. We find that the accuracy of estimation varies from stream to stream. In 
all the cases, the distributions indicate the right shape of the halo (oblate, spherical or prolate). For fitting streams B and E, the inner 
power slope 7 was kept fixed at 1.0 though its true value is 1.28. This does not introduce any significant error on the estimate of the 
density flattening. For fitting stream H, we set 7 to be a free parameter and this does not affect the estimate on q p either. For stream I, 
there are two peaks in the distribution, the actual value of q p being at the smaller peak. However, the streams generated at the smaller 
peak resemble the test stream in finer details, such as its bifurcations, which were not included in the fitting data. 



in our own galaxy with the TRGB method. In Case D, the 
distance to the progenitor and line of sight (l.o.s) velocities 
are provided in addition to the projected positions. Case E 
has the maximum information, where the rotational velocity 
curve is also available along with all the information of Case 
D. This case has been included primarily to study the effect 
of circular velocity information on the estimation of q p . The 
distributions shown have been drawn from 500,000 steps of 
the coldest Markov chain. 

Here, only the estimates of q p and V50 are presented. 
For the stream in Figure [IT] with only one prominent turn- 
ing point, the distributions of q p are limited to the oblate 
region and peak at approximately the true value of q p (red 
lines) . This particular example of this class of streams sug- 
gests that the estimates of the flattening are much more 
accurate when kinematic data are available as in cases B, 
D and E, although the estimate is good even with only the 
projection of the stream. The rightmost panels show the 
distribution of V50 for cases D and E (the cyan lines mark 
the input value). Though not shown, it is not possible to 



constrain this quantity without any kinematic information, 
as was seen for the longer streams. However, the technique 
yields very good constraints on V50 when a few velocity data 
points or the inner rotational curve are added to the spatial 
projection of the stream, making it possible to estimate the 
mass and shape of the halo for systems where such informa- 
tion is available. 

The short stream shown in Figure [18] has no notice- 
able turning point. Consequently, the distributions of q p are 
broad and highly degenerate even when the distance to the 
progenitor or line of sight velocities are provided. It is only 
when both the progenitor distance and line of sight velocities 
along the stream are provided that the distribution peaks at 
approximately the right input value of the flattening. When 
the circular velocities are also added to the fitting, then the 
distribution looks more Gaussian-like and a better estimate 
of q p is obtained. One has to bear in mind that the differ- 
ence in the estimated and true values of the flattening in 
the potential is expected to be much less than that of the 
flattening in the density distribution, since the distribution 
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Figure 12. Estimating the mass of the halo with only projections of streams. We use the circular velocity at 50kpc, V50, as an indirect 
mass parameter. The value of V50 is calculated using the potential parameters at each step of the coldest MCMC for the fitting of the 
streams shown in Figure [To] These are the distributions in V50 thus obtained for these streams. The cyan lines indicate the input V50 for 
each. It is seen that it is impossible to estimate V50 and the mass of the halo with only positional information. The algorithm finds highly 
likely solutions at different values of mass of the halo. This is because several streams (like orbits) exist that are identical in projection, 
but different in velocity space, for different values of the halo mass. As a result, the algorithm, in many cases, does not even reach the 
correct input V50 value. 



of the corresponding is narrower as per equation[9] As for 
the mass estimates, even with a stream as short as this, the 
constraints on V50 are good when sufficient information (as 
in case D) is available. Such short streams are mostly ob- 
served in systems very close to us, for which the rotational 
velocity curve, line of sight velocities and distance to the 
progenitor are available or can be easily measured. There- 
fore, they can be used as effective probes of the shapes and 
masses of nearby halos. 



6 DISCUSSION 

We have shown that stellar streams found in halos of galax- 
ies can be used to constrain parameters of the mass distri- 
bution, even without any kinematic information. The easi- 
est potential parameter to estimate, and perhaps one of the 
most interesting, is the flattening of the distribution q p . This 
is under the simplifying assumption that the halo is axisym- 
metric. In a future contribution, we will test the technique 
with triaxial halo models. Not surprisingly, our technique 



cannot constrain the mass of the halo (measured by the cir- 
cular velocity at 50 kpc) without any kinematic information. 
However, by adding the rotational velocity curve or line of 
sight velocities of the stream stars, it is possible to estimate 
the total mass as well. According to the currently-held view 
of hierarchical galaxy formation, large galaxies are built up 
by the accretion (and disruption) of smaller satellite galax- 
ies (White & Rees 19781. If this is the case, one can ex- 



pect to find remnants of the process in many galaxies and 
these will certainly be revealed by future deep surveys (E- 
ELT, JWST). By applying this technique to a large sam- 
ple of galaxies, it would be possible to obtain statistics on 
the shapes of halos, which is crucial in understanding dark 
matter and galactic evolution. The method is immediately 
applicable to individual streams, the results of which have 
several implications. For one, they would provide a better 
description of the environment in which other phenomena 
occur. The properties of the halo can be used in studies in- 
volving satellite dynamics, tidal debris, gas distributions in 
the outer regions of galaxies, lensing and many other areas 
of interest. 
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Figure 13. Estimating the outer power law f} with only projections of streams. Shown here are the distributions of obtained by fitting 
the projection of streams shown in Figure [To] The violet lines indicate the input /3 for each. This shows that /3 cannot be constrained by 
the 2-dimensional positional information alone. 




Figure 14. Effect on the estimates of parameters of adding the rotational velocity curve of the host galaxy to the projection of a stream. 
The panels from left to right show the distributions in the density flattening q p . circular velocity at 50kpc V50, and outer power law f) 
respectively, for stream B, when the inner circular velocity curve (which extends up to 30kpc for this case) is also provided in addition 
to the projected positions. The true values of each of these are marked in red, cyan and violet respectively. It is seen that the provision 
of the rotational velocity curve accurately constrains the V50 value as expected. This, in turn, helps constrain other parameters of the 
model, such as f3. 
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Figure 15. Effects on the estimates of q p and V50 when l.o.s velocities are added to the projection of the stream. The top panels show 
the distributions in q p and the bottom panels the corresponding distributions in V50 for stream G, that was generated in a perfectly 
spherical potential. The true values of q p and V50 of the model are marked in red and cyan respectively. In (a) only the projection of the 
stream on the sky is used, and a reasonably good estimate of the flattening obtained, whereas the estimate of V50 is poor. Adding the 
line of sight velocity v yo of the progenitor to the projection of the stream (but no rotational curve) greatly improves the estimates of 
V50 as well as q p (b). Adding to this five l.o.s velocity data points for each tail further improves the estimation of these parameters (c). 



The method could also provide a consistency check for 
galaxies where the dark matter distribution or its shape has 
been measured by one of the techniques mentioned in fjl] 
Tidal streams may also be used in testing alternative the- 
ories of gravity. For instance, if a prolate halo is found, it 



could eliminate MOND as a viable theory (Read & Moore 
|2005[ ). In addition to the halo parameters, the technique also 
recovers the progenitor orbit, giving us a coarse estimate of 
the distance to the progenitor, and clues to its location if 
it is unknown, as well as distances and velocities along the 
stream. 



The correction mechanism employed is one of the key 
aspects of the technique. It enables us to sample millions 
of models of the stream without having to resort to N-body 
simulations for each case, making the problem computation- 
ally tractable. During the developmental stages of this pa- 
per, we attempted to fit the tidal tails with orbits alone and 
this sometimes resulted in wrong estimates or highly degen- 
erate distributions of the flattening. Thus, the correction 
mechanism plays an important role in the effectiveness of 
the method and the accuracy of the estimates obtained. For 
its simplicity, it reproduces the stream for an orbit remark- 
ably well, even recreating the bifurcations and other features 
seen in the stream simulations. The only drawback is that 
one requires the mass of the in-falling satellite to calculate 
the corrected points, an approximate value of which can be 



obtained from star counts or integrating the light along the 
stream. 

Our correction mechanism has only been tested with 
spherical, non-rotating King satellites. It is unclear how well 
it will work for disk-like or other types of satellites, but 
it should be possible to obtain a correction mechanism for 
these based on the same principles. However, this means 
that the results for an observed stream would depend on 
the assumptions made about its progenitor. For low mass 
streams such as the ones discussed in this paper, it is safe 
to assume that the progenitor is a globular cluster or dwarf 
spheroidal galaxy. Nevertheless, it would be useful to check 
the correction mechanism with other models of the satellite 
to see how sensitive it is to the input satellite properties. 

In light of all these assumptions, the technique in its 
present form is applicable to cold, low-mass streams (M ^ 
1.0 x 10 8 Mq), for which dynamical friction can be neglected. 
Observationally, stellar streams fall into two major cate- 
gories: the ones that are closer to us, in which the positions 
and velocities of individual stars can be measured and the 
ones that are farther away and detected in surface bright- 
ness. For the latter, the kinematics are usually not mea- 
surable. In rare cases, globular clusters associated with the 
stream may be observed and their velocities can be used 
as additional information ( Mackey et al.pOlO l. Among the 
streams that are closer to us, are short structures like the 
Palomar 5 stream. For these, a large number of radial ve- 
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Figure 16. Distributions of obtained when the test streams (generated in a multicomponent host potential with a double spheroidal 
halo) are fit in a logarithmic potential, without any stellar components. The green lines indicate the flattening in potential <ju corre- 
sponding to the input flattening in density q p of the host halo calculated using equation [9] It is seen that by using a logarithmic halo to 
fit streams that were generated in a double power law density model, the flattening in the potential can still be recovered. 



locity measurements can be made, and these streams can 
be used to constrain the parameters of the host halo, de- 
spite having no turning points. Another plausible scenario is 
where two or more streams are observed in the same galaxy. 
Fitting them simultaneously would provide constraints on 
the halo, even if the individual streams are relatively short. 
Multiple streams in a halo provide a possible means of prob- 
ing the halo shape through different cross sections of the 
halo, which would also enable us to check for triaxiality of 
the halo. 

A surprising result we found, that can serve as a pow- 
erful diagnostic tool in the application of the method to real 
systems, is that if the streams (generated in a host potential 
with a double power law halo with a bulge and a disk) are fit 
in an axisymmetric logarithmic potential, the flattening of 
the potential can still be well estimated. This is encouraging 
as it suggests that while extending the method for triaxial 
models, one could use a triaxial logarithmic potential for 
preliminary estimates on the axial Battenings, which is much 
easier to use than triaxial density distributions. In addition 
to the logarithmic and double power law models described in 
this paper, we also experimented with non-parametric den- 
sity models for the halo, in which the matter distribution 
was defined by the values of densities at different radii. Pre- 



liminary tests showed that these parameters are difficult to 
constrain, but this may yet be improved upon. 



7 CONCLUSIONS 

We have extended the idea of using tidal streams to estimate 
the shape of the dark halos they reside in. The major chal- 
lenge was that for most streams, only their projected posi- 
tions are available. This has been a hindrance in using these 
streams effectively as there may be streams with similar pro- 
jections for different profiles. We overcome this difficulty by 
adopting a statistical approach. We sample various orbits 
with different orbital and potential parameters using paral- 
lel MCMC chains. The output distribution of the flattening 
parameter that we obtain from these peak at the right value, 
if the information is sufficient. We first tested the method on 
orbits in a logarithmic potential and then on streams made 
with N-body simulations in more realistic galactic potentials 
(sum of disk, ISM, bulge and spheroidal halo). Another ma- 
jor aspect by which the technique differs from earlier work 
is that the stream in question is not treated as an orbit, but 
derived as a correction to the orbit of its progenitor. 

We have demonstrated that there are cases in which it 
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Figure 17. Estimation of q p and V50 for a short stream with only one prominent turning point. The red line shows the input value of 
q p . The distributions are drawn from 500,000 steps of the coldest MCMC chain. Case A: only projected positions, red dots show the 
points used in the fitting. Case B: projected positions and l.o.s velocities at cyan squares. Case C: projected positions and distance to 
the progenitor. Case D: projected positions, distance to the progenitor and l.o.s velocities at cyan squares. Case E: Same as Case D but 
with the rotational velocity curve given. The rightmost panels show the estimation of the circular velocity at 50 kpc for case D and E, 
the cyan lines indicating its true value. It is not possible to constrain V50 without any velocity information, but if l.o.s velocities are 
provided (case D), it can be estimated even with a short stream like the one above. It is not surprising that V50 is very well constrained 
in case E as the circular velocities up to 30 kpc are given. 
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is possible to get constraints on the shape of the dark halo 
using only the projection of a tidal stream on the sky, the 
accuracy of which varies from stream to stream. The method 
is sensitive to the number of turning points of a stream and 
only works when there are at least two turning points. With 
more turning points, the accuracy of the estimate increases. 
If the stream is only a parabola or turns out to be degen- 
erate even with a few turning points, the estimate may be 
improved by adding more information such as the line of 
sight velocities or the distances along the stream if they are 
available. In our tests with simulated streams, we did not 
find cases where the fitting procedure yielded significantly 
incorrect estimates of the input halo flattening (though of 
course with insufficient information there are degenerate so- 
lutions). Exhaustive testing is required to confirm that this 
is always the case. In all our tests, we used low mass satel- 
lites {M ^ 1.0 x 10 8 M©). This method cannot be applied in 
its current form to heavy and highly disrupted streams such 
as the Giant Stellar stream around M31 ( Ibata et al.|2001a ). 
However, it may be possible to further extend the technique 
to such massive satellites by accounting for dynamical fric- 
tion as an additional correction. 

We find that the line of sight distance of the progeni- 
tor from the center of the host galaxy is well-estimated (at 
times, with a sign discrepancy) when only the projection of 
the stream on the sky is available. The other orbital param- 
eters, i.e. the velocity components of the progenitor are also 
well-estimated, but only when kinematic information, in the 
form of line of sight velocities along the stream or the H I 
rotational curve, are provided. 

We also fit the streams generated in a spheroidal den- 
sity model of the halo with orbits in a logarithmic potential. 
We see that despite using a simpler model, we are able to 
get a good estimate on the flattening of the halo. In the tests 
with the orbits in a logarithmic potential, we see that the 
initial distance yo is also constrained very well, but not the 
circular velocity. This is as expected; even for a circular or- 
bit in a Keplerian potential, it is not possible to deduce the 
central mass from the radius of the orbit alone. The method 
cannot be used to constrain the inner slope of a spheroidal 
halo as the considered halo streams are too distant to probe 
the inner mass distribution using only the projection of the 
stream. It is encouraging to note that with circular velocities 
and the projection of the stream, most of the parameters of 
the halo can be recovered, as such information is available or 
easily acquirable for several systems. This technique is read- 
ily applicable to many of the streams observed in the near 
universe such as the low mass streams around M31, NGC 
5907, NGC 891 and the Sagittarius stream in our galaxy. 
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Table 2. Models used in simulations: The parameters of the halo except its flattening, and the baryonic components of the host galaxy 
in which the streams are generated, along with their units are given, po, ro, q p , 7, /3 and rt are the central density, scale radius, density 
flattening, inner power slope, outer power slope and truncation radius of the spheroidal components (the bulge and the halo). The 
parameters of the disk components are the central surface density E d , the scalelength R d , the scalehcight z d and the central depression 
parameter R m (which is for the thin and thick disks). 



Component Parameter Model 1 Model 2 

Streams A, D, G, C, F, I Streams B, E, H 



Halo 


Po (M pc- 3 ) 
r (kpc) 

7 
P 

n (kpc) 


0.11013117 

5.23553 

1.0 

3.0 

1000 


2.7986274 

1.0 

1.28 

2.34 

1000 


Bulge 


po (M pc- 3 ) 
r (kpc) 
<?P 

7 
P 

n (kpc) 


0.65074246 

1.0 

0.6 

1.8 

1.8 

1.9 


0.3 
1.0 
0.6 
1.8 
1.8 
1.9 


Thin Disk 


S d (M oP c- 2 ) 
R d (kpc) 
z d (kpc) 


412.0 

3.2 

0.18 


341.0 

3.2 

0.18 


ISM 


E d (M Q pc- 2 ) 
Rd (kpc) 
z d (kpc) 
R m (kpc) 


69.5 
6.4 
0.04 
4.0 


57.5 
6.4 
0.04 
4.0 


Thick Disk 


E d (M pc- 2 ) 
Rd (kpc) 
z d (kpc) 


29.5 

3.2 

1.0 


24.4 

3.2 

1.0 



Table 3. Test Streams: Orbital Parameters, Halo Flattening and Age. xo,yo,zo are the initial positions and v xo ,v yo ,v zo the initial 
velocities of the orbits of the King model satellites of masses m aa t, tidal radii rt and central potential $(0)/<r 2 . The density flattening 
of the halo q p in each case is given in the first column. The ages given are the time since the starting of the simulation which gives the 
configurations of the streams that are used in the fitting. 



Stream 




xo 


yo 


zo 


v xo 






Age 


rrisat 


rt 


$(0)/<T 2 






(kpc) 


(kpc) 


(kpc) 


( km s~ 1 ) 


( kms" 1 ) 


(kms ) 


(Gyr) 


(xlO 7 M ) 


(kpc) 




A 


0.8 


30.0 


15.0 


38.0 


-110.0 


60.0 


-90.0 


7.0 


0.5 


0.8 


4.0 


B 


0.5 


35.0 


5.0 


35.0 


-100.0 


5.0 


80.0 


8.0 


0.5 


0.8 


4.0 


C 


1.6 


35.0 


20.0 


30.0 


-90.0 


80.0 


140.0 


9.0 


0.5 


0.8 


4.0 


D 


1.0 


35.0 


10.0 


35.0 


-120.0 


30.0 


100.0 


10.0 


5.0 


1.2 


4.0 


E 


0.3 


35.0 


5.0 


35.0 


-100.0 


5.0 


80.0 


5.0 


5.0 


1.5 


4.0 


F 


1.8 


35.0 


10.0 


35.0 


-120.0 


30.0 


100.0 


4.0 


5.0 


1.2 


4.0 


G 


1.0 


25.0 


20.0 


20.0 


-120.0 


90.0 


140.0 


8.0 


10.0 


1.7 


7.0 


H 


0.5 


35.0 


5.0 


35.0 


-100.0 


5.0 


80.0 


8.0 


10.0 


1.1 


2.75 


I 


2.0 


30.0 


0.0 


0.0 


180.0 


160.0 


130.0 


7.0 


10.0 


0.9 


3.0 
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